<!DOCTYPE html>
<html>
<head>
<meta content="text/html;charset=utf-8" http-equiv="Content-Type">
<meta content="utf-8" http-equiv="encoding">

<title>SPH - WebGL&WebCL Interop</title>
<meta HTTP-EQUIV="CACHE-CONTROL" CONTENT="NO-CACHE">

<style type="text/css">

#my-gui-container {
    /*margin-top: 600px;*/
}

h2 {
    color:#FFFFFF;
    margin-left: 20px;
    text-align: left;
}

p {
    color: #00FF00;
    margin-left: 20px;
}

#stats {
    position: absolute;
    top: 420px;
    width: 240px;

    /*border-style: inset;*/
    border-style: groove;
}

#submit {
    position: absolute;
    top: 520px;
    left: 10px;
    width: 240px;

    /*border-style: inset;*/
    border-style: groove;
}

#output {
    position: absolute;
    top: 0px;
    width: 0px;
    height: 0px;
    left: 720px;
    background-color: #FFF;
    display: none;
    /*border-style: groove;   */
} 

</style>

<script id="particle-vshader" type="x-shader/x-vertex">
    precision mediump float;
    uniform mat4 mvp;
    attribute vec4 position;
    attribute vec4 velocity;
    attribute vec4 acceleration;
    attribute vec2 particleIndex;
    attribute vec4 prevPosition;
    attribute vec4 sortedVelocity;
    attribute float gridCellIndex;
    attribute float gridCellIndexFixedUp;
    
    varying vec4 vColor;

    varying vec3 vVel;
    varying vec3 vPos;

    void main(void) {
        vec4 pos;
        
        pos.xyz = position.xyz;
        pos.w  = 1.0;

        vVel = velocity.xyz;

        vPos = position.xyz;

        gl_Position = mvp * pos;

        vColor = vec4(0.8, 0.8, .8, 1.0);

        gl_PointSize  = 5.0;
    }
</script>

<script id="particle-fshader" type="x-shader/x-fragment">
    
    precision mediump float;
    
    varying vec4 vColor;
    
    varying vec3 vVel;
    varying vec3 vPos;

    void main(void) {
        vec4 InnerColor = vec4(.6, .2, .10, 1.0);
        // vec4 InnerColor = vec4(0.0, 0.392157, 0.0, 1.0);
        // vec4 InnerColor = vec4(0.294118, 0.0, 0.509804, 1.0);
        vec4 OuterColor = vec4(0.0, 0.0, 0.0, 1.0);

        float dx = (gl_PointCoord.x - 0.5);
        float dy = (gl_PointCoord.y - 0.5);
        
        float r = sqrt(dx*dx + dy*dy);
        float r1 = 0.1;

        gl_FragColor = vec4(mix(InnerColor, OuterColor, smoothstep(r1, 1.0, r)));
    }

</script>

<script id="triangle-vshader" type="x-shader/x-vertex">
    precision mediump float;
    
    attribute vec3 position;
    attribute vec3 normal;

    uniform mat4 view, proj;

    varying vec3 position_eye_, normal_eye_;

    void main(void) {

        // position_eye = vec3(view * vec4(vertex_position, 1.0));
        // normal_eye = vec3(view * vec4(normal, 0.0));

        position_eye_ = vec3(view * vec4(position.xyz, 1.0));
        // normal_eye_ = vec3(view * vec4(normal, 0.0));
        // normal_eye_ = normalize(vec3(view * vec4(normal, 0.0)));
        normal_eye_ = vec3(view * vec4(normal, 0.0));
        
        gl_Position = proj * view * vec4(position.xyz, 1.0);
        
    }

</script>

<script id="triangle-fshader" type="x-shader/x-fragment">
    precision mediump float;
    
    varying vec3 position_eye_, normal_eye_;

    uniform mat4 view;

    vec3 surfaceColor = vec3(0.8, 0.1, 0.1);

    // fixed point light properties
    vec3 light_position_world  = vec3 (10.0, 5.0, 10.0);
    vec3 Ls = vec3 (0.2, 0.2, 0.2);
    vec3 Ld = vec3 (0.9, 0.9, 0.9);
    vec3 La = vec3 (1.0, 1.0, 1.0);

    // surface reflectance
    vec3 Ks = vec3 (0.5, 0.5, 0.5); // fully reflect specular light
    vec3 Kd = vec3 (0.5, 0.5, 0.5); // orange diffuse surface reflectance
    vec3 Ka = vec3 (0.1, 0.1, 0.1); // fully reflect ambient light
    float specular_exponent = 100.0; // specular 'power'

    void main(void) {
    
        // ambient intensity
        vec3 Ia = La * Ka;

        // vec3 normal = vec3(normal_eye.x+.5, normal_eye.y+.5, normal_eye.z+.5);

        // diffuse intensity
        // raise light position to eye space
        vec3 light_position_eye = vec3 (view * vec4(light_position_world, 1.0));
        vec3 distance_to_light_eye = light_position_eye - position_eye_;
        vec3 direction_to_light_eye = normalize (distance_to_light_eye);
        float dot_prod = dot (direction_to_light_eye, normalize(normal_eye_));
        dot_prod = clamp(dot_prod, 0.0, 1.0);
        vec3 Id = Ld * Kd * dot_prod * surfaceColor; // final diffuse intensity
        
        // specular intensity
        vec3 surface_to_viewer_eye = normalize(-position_eye_);
        
        // vec3 reflection_eye = reflect (-direction_to_light_eye, normalize(normal_eye));
        // float dot_prod_specular = dot (reflection_eye, surface_to_viewer_eye);
        // dot_prod_specular = max (dot_prod_specular, 0.0);
        // float specular_factor = pow (dot_prod_specular, specular_exponent);
        
        // blinn
        vec3 half_way_eye = normalize (surface_to_viewer_eye + direction_to_light_eye);
        float dot_prod_specular = clamp(dot(half_way_eye, normalize(normal_eye_)), 0.0, 1.0);
        float specular_factor = pow(dot_prod_specular, specular_exponent);
        
        vec3 Is = Ls * Ks * specular_factor * dot_prod; // final specular intensity
        
        // final colour
        gl_FragColor = vec4 (Is + Id + Ia, 1.0);
        // gl_FragColor = vec4(normalize(normal_eye_), 1.0);
        // gl_FragColor = vec4(normalize(normal_eye_), 1.0);

    }

</script>

<script id="cube-vshader" type="x-shader/x-vertex">
    precision mediump float;
    
    uniform mat4 view, proj;

    attribute vec3 cube;
    attribute vec3 vNormal;

    varying vec4 position;
    varying vec3 normal;

    void main()
    {
        position = vec3(view * vec4(cube.xyz, 1.0));

        normal = vec3(view * vec4(vNormal, 0.0));

        gl_Position = proj * view * vec4(cube.xyz, 1.0);
    }

</script>

<script id="cube-fshader" type="x-shader/x-fragment">
    precision mediump float;

    varying vec4 position;
    varying vec3 normal;

    uniform mat4 view;

    vec3 surfaceColor = vec3(0.8, 0.1, 0.1);

    // fixed point light properties
    vec3 light_position_world  = vec3 (10.0, 5.0, 10.0);
    vec3 Ls = vec3 (0.2, 0.2, 0.2);
    vec3 Ld = vec3 (0.9, 0.9, 0.9);
    vec3 La = vec3 (1.0, 1.0, 1.0);

    // surface reflectance
    vec3 Ks = vec3 (0.5, 0.5, 0.5); // fully reflect specular light
    vec3 Kd = vec3 (0.5, 0.5, 0.5); // orange diffuse surface reflectance
    vec3 Ka = vec3 (0.1, 0.1, 0.1); // fully reflect ambient light
    float specular_exponent = 100.0; // specular 'power'

    void main() {
        
        // ambient intensity
        vec3 Ia = La * Ka;

        // vec3 normal = vec3(normal_eye.x+.5, normal_eye.y+.5, normal_eye.z+.5);

        // diffuse intensity
        // raise light position to eye space
        vec3 light_position_eye = vec3 (view * vec4(light_position_world, 1.0));
        vec3 distance_to_light_eye = light_position_eye - position;
        vec3 direction_to_light_eye = normalize (distance_to_light_eye);
        float dot_prod = dot (direction_to_light_eye, normalize(normal));
        dot_prod = clamp(dot_prod, 0.0, 1.0);
        vec3 Id = Ld * Kd * dot_prod * surfaceColor; // final diffuse intensity
        
        // specular intensity
        vec3 surface_to_viewer_eye = normalize(-position);
        
        // vec3 reflection_eye = reflect (-direction_to_light_eye, normalize(normal_eye));
        // float dot_prod_specular = dot (reflection_eye, surface_to_viewer_eye);
        // dot_prod_specular = max (dot_prod_specular, 0.0);
        // float specular_factor = pow (dot_prod_specular, specular_exponent);
        
        // blinn
        vec3 half_way_eye = normalize (surface_to_viewer_eye + direction_to_light_eye);
        float dot_prod_specular = clamp(dot(half_way_eye, normalize(normal)), 0.0, 1.0);
        float specular_factor = pow(dot_prod_specular, specular_exponent);
        
        vec3 Is = Ls * Ks * specular_factor * dot_prod; // final specular intensity
        
        // final colour
        // gl_FragColor = vec4 (Is + Id + Ia, 1.0);
        gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0);
    }
</script>

<script id="plane-vshader" type="x-shader/x-vertex">

    precision mediump float;
    uniform mat4 mvp;
    uniform mat4 np;

    attribute vec3 plane;
    attribute vec3 vNormal;

    varying vec4 position;
    varying vec3 light;
    varying vec3 normal;
    
    void main()
    {
        position = vec4(plane.xyz, 1.0);
        position = mvp * position;
        
        gl_Position = position;
    }   

</script>

<script id="plane-fshader" type="x-shader/x-fragment">

    precision mediump float;

    varying vec4 position;
    varying vec3 light;
    varying vec3 normal;

    void main() {
        gl_FragColor = vec4(0.5, 0.08, 0.1, 1.0);

    }

</script>

<script id="sph_kernel_applyBodyForce" type="x-kernel">

    __kernel void sph_kernel_applyBodyForce(
            int kParticleCount,
            __global float4* vel,
            float dt)
    {
        unsigned int gid = get_global_id(0);

        // vel[gid].y -= 9.8 * dt;
        vel[gid].y -= 13.546 * dt;
    }

</script>

<script id="sph_kernel_advance" type="x-kernel">

    __kernel void sph_kernel_advance(
        int kParticleCount,
        float dt,
        __global float4* position,
        __global float4* prevPos,
        __global float4* velocity)
    {
        unsigned int gid = get_global_id(0);
        
        prevPos[gid].xyzw = position[gid].xyzw;
        position[gid].xyz += (float3)(dt) * velocity[gid].xyz;
    }

</script>

<script id="sph_kernel_hashparticles" type="x-kernel">
    
    __kernel void sph_kernel_hashparticles(
        float viewWidth,
        float viewHeight,
        float viewDepth,
        int nx,
        int ny,
        int nz,
        float cellSize,
        __global float4* pos,
        __global uint2* partIdx,
        __global float4* vel)
    {
        unsigned int gid = get_global_id(0);

        float3 tempPos;
        // pos(x, y, z) in [-viewWidth/2, viewWidth/2] mapped to [0, nx*cellSize]
        // newPos = (val - src0) / (src1 - src0) * (dst1 - dst0) + dst0
        float rangeRatioX = (nx * cellSize - 0.0f) / viewWidth;
        float rangeRatioY = (ny * cellSize - 0.0f) / viewHeight;
        float rangeRatioZ = (nz * cellSize - 0.0f) / viewDepth;
        
        tempPos.x = (pos[gid].x - (-viewWidth/2) ) * rangeRatioX;
        tempPos.y = (pos[gid].y - (-viewHeight/2)) * rangeRatioY;
        tempPos.z = (pos[gid].z - (-viewDepth/2) ) * rangeRatioZ;

        // myPos(x, y, z) in [0, nx*cellSize] mapped to 3d voxel coords [v3d.x, v3d.y, v3d.z]
        // each between the range [0, n] (assuming cubic space)
        float3 voxel3d;
        voxel3d.x = tempPos.x / cellSize;
        voxel3d.y = tempPos.y / cellSize;
        voxel3d.z = tempPos.z / cellSize;

        int voxelID = floor(voxel3d.x) + floor(voxel3d.y) * nx + floor(voxel3d.z) * nx * ny;
        
        pos[gid].w = voxelID;
        
        partIdx[gid].x = voxelID;
        partIdx[gid].y = gid;        

    }

</script>

<!-- 
    Bitonic sort code source: 
    https://github.com/clockfort/amd-app-sdk-fixes/blob/master/samples/opencl/cl/app/BitonicSort/BitonicSort_Kernels.cl 
-->
<script id="sph_kernel_sort" type="x-kernel">

    __kernel void sph_kernel_sort (
        __global uint2* partIdx,
        int stage,
        int passOfStage,
        int direction) 
    {
        uint sortIncreasing = direction;
        uint threadId = get_global_id(0);
    
        uint pairDistance = 1 << (stage - passOfStage);
        uint blockWidth   = 2 * pairDistance;

        uint leftId = (threadId % pairDistance) 
                   + (threadId / pairDistance) * blockWidth;

        uint rightId = leftId + pairDistance;
    
        uint2 leftElement = partIdx[leftId];
        uint2 rightElement = partIdx[rightId];
    
        uint sameDirectionBlockWidth = 1 << stage;
    
        if((threadId/sameDirectionBlockWidth) % 2 == 1)
            sortIncreasing = 1 - sortIncreasing;

        uint2 greater;
        uint2 lesser;
        if(leftElement.x > rightElement.x)
        {
            greater = leftElement;
            lesser  = rightElement;
        }
        else
        {
            greater = rightElement;
            lesser  = leftElement;
        }
    
        if(sortIncreasing)
        {
            partIdx[leftId]  = lesser;
            partIdx[rightId] = greater;
        }
        else
        {
            partIdx[leftId]  = greater;
            partIdx[rightId] = lesser;
        }
    }

</script>

<script id="sph_kernel_sortPostPass" type="x-kernel">
    
    __kernel void sph_kernel_sortPostPass (
        __global float4* pos,
        __global float4* vel,
        __global uint2* partIdx,
        __global float4* sortedPos,
        __global float4* sortedVel,
        __global float4* prevPos,
        __global float4* sortedPrevPos)
    {
        unsigned int gid = get_global_id(0);

        int particleID = partIdx[gid].y;

        // float4 tempPos = pos[particleID];

        sortedPos[gid].x = pos[particleID].x;
        sortedPos[gid].y = pos[particleID].y;
        sortedPos[gid].z = pos[particleID].z;
        sortedPos[gid].w = pos[particleID].w;

        // float4 tempVel = vel[particleID];

        sortedVel[gid].x = vel[particleID].x;
        sortedVel[gid].y = vel[particleID].y;
        sortedVel[gid].z = vel[particleID].z;
        sortedVel[gid].w = vel[particleID].w;

        // float4 tempPrevPos = prevPos[particleID];

        sortedPrevPos[gid].x = prevPos[particleID].x;
        sortedPrevPos[gid].y = prevPos[particleID].y;
        sortedPrevPos[gid].z = prevPos[particleID].z;
        sortedPrevPos[gid].w = prevPos[particleID].w;

        barrier(CLK_GLOBAL_MEM_FENCE);

        // pos[gid].x = sortedPos[gid].x;
        // pos[gid].y = sortedPos[gid].y;
        // pos[gid].z = sortedPos[gid].z;
        // pos[gid].w = sortedPos[gid].w;

        // vel[gid].x = sortedVel[gid].x;
        // vel[gid].y = sortedVel[gid].y;
        // vel[gid].z = sortedVel[gid].z;
        // vel[gid].w = sortedVel[gid].w;
        
        pos[particleID].x = sortedPos[particleID].x;
        pos[particleID].y = sortedPos[particleID].y;
        pos[particleID].z = sortedPos[particleID].z;
        pos[particleID].w = sortedPos[particleID].w;

        barrier(CLK_GLOBAL_MEM_FENCE);

        vel[particleID].x = sortedVel[particleID].x;
        vel[particleID].y = sortedVel[particleID].y;
        vel[particleID].z = sortedVel[particleID].z;
        vel[particleID].w = sortedVel[particleID].w;
    }

</script>

<script id="sph_kernel_indexx" type="x-kernel">

    __kernel void sph_kernel_indexx(
        int kParticleCount,
        // __global float4* sortedPos,
        __global float4* pos,
        __global int* gridCellIdx)
    {
        unsigned int gid = get_global_id(0);
        gridCellIdx[gid] = -1;

        // binary search into sortedPos and find first particle
        // in voxel gid

        int low = 0;
        int hi = kParticleCount - 1;
        int mid = 0;

        while (low <= hi) {
            mid = (hi + low) / 2;
            // if (floor(sortedPos[mid].w) == gid) {
            if (floor(pos[mid].w) == gid) {    
                int front = mid - 1;
                
                // while((front >= 0) && (floor(sortedPos[front].w) == gid))
                while((front >= 0) && (floor(pos[front].w) == gid))
                {
                    front--;
                }

                if (mid > front - 1)
                {
                    gridCellIdx[gid] = front + 1;
                }

                break;

            // } else if ( sortedPos[mid].w < gid) {
            } else if ( pos[mid].w < gid) {                
                low = mid + 1;
            } else {
                hi = mid - 1;
            }
        }
    }

</script>

<script id="sph_kernel_indexPostPass" type="x-kernel">
    
    __kernel void sph_kernel_indexPostPass(
        int kParticleCount,
        __global int* gridCellIdx,
        __global int* gridCellIdxFixedUp)
    {
        unsigned int gid = get_global_id(0);
        
        if (gridCellIdx[gid] != -1)
        {
            gridCellIdxFixedUp[gid] = gridCellIdx[gid];
        } else {
            int preCell = gid;
        
            while (preCell >= 0) 
            {
                int pid = gridCellIdx[preCell];
                if (pid != -1)
                {
                    gridCellIdxFixedUp[gid] = pid;
                    break;
                }
                else {
                    preCell--;
                }
            }
        }
    }

</script>

<script id="sph_kernel_findNeighbors" type="x-kernel">

    __kernel void sph_kernel_findNeighbors(
        int kParticleCount,
        float viewW,
        float viewH,
        float viewD,
        int nx,
        int ny,
        int nz,
        float cellSize,
        __global float4* pos,
        __global int* gridCellIdxFixedUp,
        __global int* neighborMap,
        int maxNeighbors,
        float kEpsilon)
    {
        unsigned int gid = get_global_id(0);
        
        float3 myPos = pos[gid].xyz;

        int count = 0;
        int i;
        int j;

        // set all entries to be the particle itself
        for (i = 0; i < maxNeighbors; i++)
        {
            neighborMap[maxNeighbors*gid+i] = gid;
        }

        for (j = 0; j < kParticleCount; j++)
        {
            // float3 particleJ = sortedPos[j].xyz;
            float3 particleJ = pos[j].xyz;
            float dist = length(myPos.xyz - particleJ);

            if ((dist >= kEpsilon) && (dist <= cellSize) && (j != gid))
            {
                neighborMap[maxNeighbors*gid+count] = j;
                count++;
                if (count == maxNeighbors)
                    break;
            }
        }
    }

    // __kernel void sph_kernel_findNeighbors(
    //     int kParticleCount,
    //     float viewW,
    //     float viewH,
    //     float viewD,
    //     int nx,
    //     int ny,
    //     int nz,
    //     float cellSize,
    //     __global float4* pos,
    //     __global int* gridCellIdxFixedUp,
    //     __global int* neighborMap,
    //     int maxNeighbors,
    //     float kEpsilon)
    // {
    //     unsigned int gid = get_global_id(0);
        
    //     float3 myPos = pos[gid].xyz;

    //     float rangeRatioX = (nx * cellSize - 0) / viewW;
    //     float rangeRatioY = (ny * cellSize - 0) / viewH;
    //     float rangeRatioZ = (nz * cellSize - 0) / viewD;
        
    //     float4 tempPos;
    //     tempPos.x = (myPos.x - (-viewW/2)) * rangeRatioX + 0;
    //     tempPos.y = (myPos.y - (-viewH/2)) * rangeRatioY + 0;
    //     tempPos.z = (myPos.z - (-viewD/2)) * rangeRatioZ + 0;

    //     // myPos(x, y, z) in [0, nx*2*h] mapped to 3d voxel coords [v3d.x, v3d.y, v3d.z]
    //     // each between the range [0, n] (assuming cubic space)
    //     float3 voxel3d;
    //     voxel3d.x = tempPos.x / cellSize;
    //     voxel3d.y = tempPos.y / cellSize;
    //     voxel3d.z = tempPos.z / cellSize;

    //     float3 corner;
        
    //     corner.x = floor(voxel3d.x + .5);
    //     corner.y = floor(voxel3d.y + .5);
    //     corner.z = floor(voxel3d.z + .5);

    //     // int voxelID = floor(voxel3d.x) + floor(voxel3d.y) * nx + floor(voxel3d.z) * nx * ny;

    //     // compute voxel IDs based on the position of the particles 
    //     int voxel[8];
    //     voxel[0] = corner.x + corner.y * nx + corner.z * nx * ny ;
    //     voxel[1] = ((corner.x < 1) ? 0 : (corner.x-1)) + corner.y * nx + corner.z * nx * ny;
    //     voxel[2] = corner.x + ((corner.y < 1) ? 0 : (corner.y - 1)) * nx + corner.z * nx * ny;
    //     voxel[3] = ((corner.x < 1) ? 0 : (corner.x-1)) + ((corner.y < 1) ? 0 : (corner.y - 1)) * nx + corner.z * nx * ny;
    //     voxel[4] = corner.x + corner.y * nx + ((corner.z < 1) ? 0 : (corner.z - 1)) * nx * ny;
    //     voxel[5] = ((corner.x < 1) ? 0 : (corner.x-1)) + corner.y * nx + ((corner.z < 1) ? 0 : (corner.z - 1)) * nx * ny;
    //     voxel[6] = corner.x + ((corner.y < 1) ? 0 : (corner.y - 1)) * nx + ((corner.z < 1) ? 0 : (corner.z - 1)) * nx * ny;
    //     voxel[7] = ((corner.x < 1) ? 0 : (corner.x-1)) + ((corner.y < 1) ? 0 : (corner.y - 1)) * nx + ((corner.z < 1) ? 0 : (corner.z - 1)) * nx * ny;

    //     int random_voxel[8] = {0, 1, 2, 3, 4, 5, 6, 7};

    //     // // use work item ID, gid, so seed uvec2 v = (gid, gid+1) etc.
    //     // uint2 rvoxel01 = uint2(gid, gid+1);
    //     // uint2 rvoxel02 = uint2(gid+2, gid+3);
    //     // uint2 rvoxel03 = uint2(gid+4, gid+5);
    //     // uint2 rvoxel04 = uint2(gid+6, gid+7);
        
    //     // rvoxel01 = TEA8(rvoxel01);
    //     // rvoxel02 = TEA8(rvoxel02);
    //     // rvoxel03 = TEA8(rvoxel03);
    //     // rvoxel04 = TEA8(rvoxel04);

    //     // float2 rvoxel01F = rand(rvoxel01);
    //     // float2 rvoxel02F = rand(rvoxel02);
    //     // float2 rvoxel03F = rand(rvoxel03);
    //     // float2 rvoxel04F = rand(rvoxel04);

    //     // // rangeMap(float val, float src0, float src1, float dst0, float dst1);
    //     // rvoxel01F.x = rangeMap(rvoxel01F.x, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel01F.y = rangeMap(rvoxel01F.y, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel02F.x = rangeMap(rvoxel02F.x, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel02F.y = rangeMap(rvoxel02F.y, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel03F.x = rangeMap(rvoxel03F.x, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel03F.y = rangeMap(rvoxel03F.y, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel04F.x = rangeMap(rvoxel04F.x, 0.0, 1.0, 0.0, 8.0);
    //     // rvoxel04F.y = rangeMap(rvoxel04F.y, 0.0, 1.0, 0.0, 8.0);
        
    //     // int random_voxel[8] = {
    //     //     (int)rvoxel01F.x, (int)rvoxel01F.y,
    //     //     (int)rvoxel02F.x, (int)rvoxel02F.y,
    //     //     (int)rvoxel03F.x, (int)rvoxel03F.y,
    //     //     (int)rvoxel04F.x, (int)rvoxel04F.y};
        
    //     int i;
    //     int j;
    //     int count = 0;
    //     // find 32 neighboring particles
    //     // not handling "not enough neighboring particle" situation
    //     // for now, assume we can find 32 neighbors for sure
    //     for (i = 0; i < 8; i++)
    //     {
    //         // locate the voxel currently being investigated
    //         int vID = voxel[random_voxel[i]];

    //         int pid = gridCellIdxFixedUp[vID];
                
    //         // in voxel vID get a random number based on the number of particles in that voxel
    //         // first, get the number of particles in voxel[voxelID], pCount
    //         int npid = (vID + 1 < nx * ny * nz) ? gridCellIdxFixedUp[vID+1] : kParticleCount;
    //         int pCount = npid - pid;

    //         // if the voxel is empty, skip this iteration
    //         if (pCount == 0)
    //         {
    //             continue;
    //         }

    //         // uint2 rparticle = uint2(pCount, pCount+1);
    //         // rparticle = TEA8(rparticle);
    //         // float2 rparticleF = rand(rparticle);
    //         // rparticleF.x = rangeMap(rparticleF.x, 0.0, 1.0, 0.0, (float)pCount);
    //         // // rparticleF.y = rangeMap(rparticleF.y, 0.0, 1.0, 0.0, (float)pCount);

    //         // int randPart = (int)rparticleF.x;
            
    //         for (j = pid; j < pid + pCount; j++)
    //         {
    //             if (count < maxNeighbors)
    //             {
    //                 float3 diff = myPos.xyz - pos[j].xyz;
    //                 float dist2 = diff.x * diff.x + diff.y * diff.y + diff.z * diff.z;
                    
    //                 if (dist2 < cellSize * cellSize)
    //                 {
    //                     // int tempIdx = j;
    //                     neighborMap[maxNeighbors*gid+count] = j;
    //                     count++;
    //                 }
    //                 else
    //                     continue;
    //             }
    //         }
    //     }
    // }

</script>

<script id="sph_kernel_pressure" type="x-kernel">

    __kernel void sph_kernel_pressure(
        float m,
        float cellSize,
        float kNorm,
        float kNearNorm,
        int maxNeighbors,
        float kStiffness,
        float kNearStiffness,
        float kRestDensity,
        float kEpsilon,
        __global float4* pos,
        __global int* neighborMap,
        __global float* density,
        __global float* nearDensity,
        __global float* pressure,
        __global float* nearPressure)
    {
        unsigned int gid = get_global_id(0);    

        float3 myPos = pos[gid].xyz;

        int i;
        float3 posJ;
        float tempDens = 0.0;
        float tempNearDens = 0.0;

        for (i = 0; i < maxNeighbors; i++)
        {
            int nID = neighborMap[maxNeighbors * gid+i];
            if (nID == gid)
                continue;

            posJ = pos[nID].xyz;

            float r = length(myPos - posJ);

            float a = 1 - r/cellSize;

            tempDens += m * a * a * a * kNorm;
            tempNearDens += m * a * a * a * a * kNearNorm;
            
        }

        density[gid] = tempDens;
        nearDensity[gid] = tempNearDens;
        pressure[gid] = kStiffness * (tempDens - m * kRestDensity);
        nearPressure[gid] = kNearStiffness * tempNearDens;

    }

</script>

<script id="sph_kernel_calcRelaxPos" type="x-kernel">

    __kernel void sph_kernel_calcRelaxPos(
        float m,
        float cellSize,
        float dt,
        float dt2,
        float kNearNorm,
        float kNorm,
        float kSurfaceTension,
        float kLinearViscocity,
        float kQuadraticViscocity,
        int maxNeighbors,
        __global float4* pos,
        __global float4* vel,
        __global int* neighborMap,
        __global float* density,
        __global float* nearDensity,
        __global float* pressure,
        __global float* nearPressure,
        __global float4* relaxPos)
    {
        unsigned int gid = get_global_id(0);

        float3 myPos = pos[gid].xyz;
        float3 myVel = vel[gid].xyz;

        float x = myPos.x;
        float y = myPos.y;
        float z = myPos.z;

        int i;
        float3 posJ;
        float3 velJ;

        for (i = 0; i < maxNeighbors; i++)
        {
            int nID = neighborMap[maxNeighbors * gid+i];
            if (nID == gid)
                continue;

            posJ = pos[nID].xyz;

            float3 diff = posJ - myPos;
            float r = length(diff);

            float dx = diff.x;
            float dy = diff.y;
            float dz = diff.z;

            float a = 1 - r/cellSize;

            float d = dt2 * ((nearPressure[gid] + nearPressure[nID]) * a * a * a * kNearNorm + (pressure[gid] + pressure[nID]) * a * a * kNorm) / 2.0;

            x -= d * dx / (r * m);
            y -= d * dy / (r * m);
            z -= d * dz / (r * m);

            x += (kSurfaceTension/m) * m * a * a * kNorm * dx;
            y += (kSurfaceTension/m) * m * a * a * kNorm * dy;
            z += (kSurfaceTension/m) * m * a * a * kNorm * dz;

            velJ = vel[nID].xyz;

            float3 diffV = myVel - velJ;
            float u = diffV.x * dx + diffV.y * dy + diffV.z * dz;

            if (u > 0)
            {
                u /= r;

                float a = 1 - r/cellSize;

                float I = .5 * dt * a * (kLinearViscocity * u + kQuadraticViscocity * u * u);

                x -= I * dx * dt;
                y -= I * dy * dt;
                z -= I * dz * dt;
            }
        }

        relaxPos[gid].x = x;
        relaxPos[gid].y = y;
        relaxPos[gid].z = z;

    }

</script>

<script id="sph_kernel_moveToRelaxPos" type="x-kernel">

    __kernel void sph_kernel_moveToRelaxPos(
        float dt,
        __global float4* pos,
        __global float4* prevPos,
        __global float4* vel,
        __global float4* relaxPos)
    {
        unsigned int gid = get_global_id(0);

        pos[gid].x = relaxPos[gid].x;
        pos[gid].y = relaxPos[gid].y;
        pos[gid].z = relaxPos[gid].z;

        vel[gid].x = (pos[gid].x - prevPos[gid].x) / dt;
        vel[gid].y = (pos[gid].y - prevPos[gid].y) / dt;
        vel[gid].z = (pos[gid].z - prevPos[gid].z) / dt;

    }

</script>

<script id="sph_kernel_resolveCollisions" type="x-kernel">

    __kernel void sph_kernel_resolveCollisions(
        float dt,
        float pRadius,
        float cellSize,
        float viewW,
        float viewH,
        float viewD,
        __global float4* pos,
        __global float4* vel)
    {
        unsigned int gid = get_global_id(0);

        float3 myPos = pos[gid].xyz;
        float3 myVel = vel[gid].xyz;

        float3 center         = float3(0.0);
        float3 boxSize        = float3(viewW/2-cellSize-pRadius, viewH/2-cellSize-pRadius, viewD/2-cellSize-pRadius);
        
        float3 xLocal = myPos - center;
        float3 contactPointLocal = min(boxSize, max(-boxSize, xLocal));
        float3 contactPoint = contactPointLocal + center;
        float distance = length(contactPoint - myPos);
        
        if (distance > 0.0 && length(myVel) > 0.0) 
        {
            float3 normal = normalize(sign(contactPointLocal - xLocal));
            float restitution = .5*distance / (dt * length(myVel));
            
            vel[gid].xyz -= (float3)((1.0 + restitution) * dot(myVel, normal)) * normal;
            pos[gid].xyz = contactPoint;
        }

    }

</script>

<script id="mc_kernel_reset" type="x-kernel">
    __kernel void mc_kernel_reset(
        int numVolIdx,
        __global float4* mcgrid)
    {
        unsigned int gid = get_global_id(0);

        mcgrid[gid].w = 0.0;
    }

</script>

<script id="mc_kernel_gridval" type="x-kernel">
    
    __kernel void mc_kernel_gridval(
        float m,
        float cellSize,
        float kNorm,
        float kNearNorm,
        float kEpsilon,
        __global float4* pos,
        __global float4* mcgrid,
        int numVolIdx)
    {
        unsigned int gid = get_global_id(0);

        float3 myPos = pos[gid].xyz;

        int i;
        float3 posJ;
        
        float a;
        
        for (i = 0; i < numVolIdx; i++)
        {
            posJ = mcgrid[i].xyz;
            a = 0.0;

            float r = length(myPos - posJ);

            if (r >= kEpsilon && r <= cellSize)
            {
                // a = 1.0 - r/cellSize;
                // a = r/cellSize;
                a = cellSize/r;
            }
            else
                continue;

            // mcgrid[i].w += m * (.5*cos(a)) * kNearNorm;
            mcgrid[i].w += m * a * kNearNorm;
            // mcgrid[i].w += m * a * kNorm;
        }
    }

</script>

<script id="mc_kernel_gridGrad" type="x-kernel">

    __kernel void mc_kernel_gridGrad(
        int volnx,
        int volny,
        int volnz,
        float edgex,
        float edgey,
        float edgez,
        __global float4* mcgrid,
        __global float4* mcgridGrad)
    {
        unsigned int gid = get_global_id(0);
        
        if (mcgridGrad[gid].w != -1.0)
        {
            mcgridGrad[gid].x = (mcgrid[gid-1].w - mcgrid[gid+1].w)/edgex;
            mcgridGrad[gid].y = (mcgrid[gid-(volny+1)].w - mcgrid[gid+(volny+1)].w)/edgey;
            mcgridGrad[gid].z = (mcgrid[gid-(volny+1)*(volnx+1)].w - mcgrid[gid+(volny+1)*(volnx+1)].w)/edgez;
        }
    }

</script>

<script id="mc_kernel_cubeindex" type="x-kernel">
    
    __kernel void mc_kernel_cubeindex(
        float isoval, // isosurface value
        int volnx,
        int volny,
        int volnz,
        __global float4* mcgrid,
        __global int* pointindex,
        __global int* cubeindex)
    {
        unsigned int gid = get_global_id(0);

        if ((int)gid < volnx * volny * volnz)
        {
            int idx0 = pointindex[8*gid+0];
            int idx1 = pointindex[8*gid+1];
            int idx2 = pointindex[8*gid+2];
            int idx3 = pointindex[8*gid+3];
            int idx4 = pointindex[8*gid+4];
            int idx5 = pointindex[8*gid+5];
            int idx6 = pointindex[8*gid+6];
            int idx7 = pointindex[8*gid+7];

            float val0 = mcgrid[idx4].w;
            float val1 = mcgrid[idx5].w;
            float val2 = mcgrid[idx1].w;
            float val3 = mcgrid[idx0].w;
            float val4 = mcgrid[idx7].w;
            float val5 = mcgrid[idx6].w;
            float val6 = mcgrid[idx2].w;
            float val7 = mcgrid[idx3].w;

            int cidx = 0;

            cidx =  (val0 < isoval); 
            cidx += (val1 < isoval)*2; 
            cidx += (val2 < isoval)*4; 
            cidx += (val3 < isoval)*8; 
            cidx += (val4 < isoval)*16; 
            cidx += (val5 < isoval)*32; 
            cidx += (val6 < isoval)*64; 
            cidx += (val7 < isoval)*128;

            cubeindex[gid] = cidx;
        }
    }
</script>

<script type="text/javascript" src="/socket.io/socket.io.js"></script> 
<!--The above library(/socket.io/socket.io.js) will be generated by socket.io module of server -->
<script type="text/javascript" src="http://code.jquery.com/jquery-1.10.2.min.js"></script>

<script type="text/javascript" src="Util/dat.gui.js"></script>
<script src="Util/sampler.js"></script>
<script src="Util/shader.js"></script>
<script src="Util/J3DIMath.js"></script>

<script src="Controller.js"></script>
<script src="InitSim.js"></script>
<script src="MarchingCubes.js"></script>
<script src="SimCL.js"></script>
<script src="SimJS.js"></script>
<script src="DrawGL.js"></script>
<script src="DrawJS.js"></script>

</head>

<body onload="onLoad()" bgColor="#202020">
    <!-- <div style="position:absolute; left:0px; top:0px">
        <div id="sim" class="info" style="position:absolute; width: 0px;"></div>
        <div id="sms" class="info" style="position:absolute; width: 0px;"></div>
        <div id="drw" class="info" style="position:absolute; width: 0px;"></div>
        <div id="dms" class="info" style="position:absolute; width: 0px;"></div>

    </div> -->
    
    <div id='my-gui-container'></div>

    <!-- canvas must be square since simulator works in normalized device coordinates -->
    <canvas id="canvas2D" style="z-index: -1; position:absolute; left:300px; top:0px;"></canvas>
    <canvas id="canvas3D" style="z-index: -1; position:absolute; left:300px; top:0px;"></canvas>
    
    <div id='stats'>
        <h2>Default settings</h2>
        <p>OpenCL device: CPU</p>
    </div>

    <div id="output"></div>
    <form>

        <!-- <input type="text" id="message" /> text form to send data to the server -->
        <!-- <input id="submit" type="button" value="Send data to Server"> -->
        <div id="content"></div> <!--This is where the data from the server is added-->
     
    </form>
    
    <script type="text/javascript">
                var socket = io.connect("/"); 
                
                $(function(){
                    $('#submit').click(function(){ /*listening to the button click using Jquery listener*/
                        
                        var data = document.getElementById("output").innerHTML;
                        socket.send(data);
     
                    });
                });
     
    </script>
    
</body>
</html>